function [M, L, info] = kmedoid_std(C, M0, varargin)
%KMEDOID_STD Implements standard K-medoid clustering
%
%   [M, L] = KMEDOID_STD(C, M0, ...);
%   [M, L] = KMEDOID_STD(C, K, ...);
%   [M, L] = KMEDOID_STD({cfun, n}, M0, ...);
%   [M, L] = KMEDOID_STD({cfun, n}, K, ...);
%
%       performs K-medoid clustering based on the pairwise cost matrix C.
%
%       Suppose there are n samples. 
%
%       Input arguments:
%       - C:        The cost matrix, of size n x n. 
%                   C(i, j) is the cost of assigning i-th sample to the
%                   cluster with the j-th sample as the center.
%
%       - cfun:     The pairwise cost function.
%
%       - n:        The number of samples
%
%       - M0:       The initial selected medoids (an index vector of 
%                   length K > 1)
%
%       - K:        The number of medoids. 
%
%       If K instead of M0 is given, then M0 will be generated by randomly
%       picking K samples. One can also use kmseeds to pick the initial
%       medoids.
%
%       Output arguments:
%       - M:        The resultant medoids (a row vector of indices)
%       - L:        The assignment of samples to medoids [1 x n].
%
%       Options:
%       
%       The user can specify options (in name/value pairs) to 
%       control the algorithm. The options are listed as below:
%
%       - MaxIter:  the maximum number of iterations (default = 100)
%
%       - TolC:     the maximum allowable number of label changes at
%                   convergence. (default = 0)
%
%       - TolFun:   the maximum allowable change of objective value at
%                   convergence (default = 1e-8)
%
%       - Display:  the level of displaying (default = 'off')
%                   - 'off':    display nothing
%                   - 'final':  display a brief summary when the procedure
%                               finishes
%                   - 'iter':   display at each iteration
%                  
%       The user can also use kmopts function to construct
%       an option struct and use it as an input argument in the place
%       of the name/value pairs.
%
%
%   [M, L, info] = KMEDOID_STD( ... );
%       
%       Additionally returns a struct with the information of the procedure
%       - niters:       The number of iterations
%       - objv:         The objective value (total costs)
%       - converged:    Whether the procedure converged.

%   History
%       - Created by Dahua Lin, on Jun 4, 2008
%       - Modified by Dahua Lin, on Sep 28, 2010
%       - Modified by Dahua Lin, on Mar 28, 2012
%

%% parse and verify input

if isnumeric(C)
    n = size(C, 1);
    if ~(isfloat(C) && isreal(C) && ndims(C) == 2 && size(C,2) == n)
        error('kmedoid_std:invalidarg', ...
            'C should be a real square matrix.');
    end
    use_mat = true;
    
elseif iscell(C) && numel(C) == 2    
    n = C{2};
    C = C{1};    
    if ~(isa(C, 'function_handle') && isnumeric(n) && isscalar(n))
        error('kmedoid:invalidarg', 'The first argument is invalid.');
    end    
    use_mat = false;
    
else
    error('kmedoid_std:invalidarg', 'The first argument is invalid.');
end

if isscalar(M0)
    K = M0;
    M0 = [];
    
elseif isnumeric(M0) && isvector(M0)
    K = numel(M0);
    if size(M0, 1) > 1
        M0 = M0.';
    end
    
else
    error('kmedoid_std:invalidarg', 'The second argument is invalid.');
end

if ~(isnumeric(K) && (K >= 1 && K <= n) && K == fix(K))
    error('kmedoid_std:invalidarg', 'The value of K is invalid.');
end

% get options

opts = kmopts(varargin{:});

tolfun = opts.tolfun;
tolc = opts.tolc;
maxiter = opts.maxiter;
displevel = opts.displevel;

%% main 

% initialize medoids

if isempty(M0)    
    M0 = randpick(n, K); 
end
M = M0;

% initialize assignments and objective

costs = calc_costs(C, M, n, use_mat);
[min_costs, L] = min(costs, [], 2);
L = L.';

objv = sum(min_costs);

% main loop

t = 0;
converged = false;

if displevel >= 2
    print_iter_header();
end


while ~converged && t < maxiter
    
    t = t + 1;
    
    % identify affected clusters
    
    if t == 1
        aff_cs = 1 : K;
    else        
        aff_cs = find(intcount(K, [L(chs), L_pre(chs)]));
    end
                      
    % update centers        
            
    gs = intgroup(K, L);
    for k = aff_cs
        gk = gs{k};
        Ck = C(gk, gk);
        [~, si] = min(sum(Ck, 1));
        M(k) = gk(si);
    end        
        
    % re-compute costs
    
    if numel(aff_cs) >= K / 2
        costs = calc_costs(C, M, n, use_mat);
    else
        costs(:, aff_cs) = calc_costs(C, M(:, aff_cs), n, use_mat);
    end
    
    % re-assign samples to centers
    L_pre = L;
    [min_costs, L] = min(costs, [], 2);
    L = L.';
    
    % determine convergence
    
    objv_pre = objv;
    objv = sum(min_costs);
    
    v_ch = objv - objv_pre;
    
    chs = find(L ~= L_pre);
    converged = (abs(v_ch) <= tolfun || numel(chs) <= tolc);
    
    % print iteration info
    
    if displevel >= 2
        print_iter(t, numel(chs), aff_cs, objv, v_ch);
    end    

end

if displevel >= 1
    print_final(t, objv, converged);
end

if nargout >= 3
    info.niters = t;
    info.objv = objv;
    info.converged = converged;
end


%% Auxiliary functions

function costs = calc_costs(C, M, n, use_mat)

if use_mat
    costs = C(:, M);
else
    costs = C(1:n, M);
end


%% Printing functions

function print_iter_header()

fprintf(' Iter    # ch.assign (aff.clus)     objv (change)  \n');
fprintf('------------------------------------------------------\n');


function print_iter(it, ch, afc, objv, vch)

fprintf('%5d        %7d (%5d)   %12.6g (%.4g)\n', it, ch, numel(afc), objv, vch);


function print_final(it, objv, converged)

fprintf('K-medoid: total_cost = %.6g, converged = %d [total # iters = %d]\n', ...
    objv, converged, it);



